package scu.maqiang.numeric;

import scu.maqiang.mesh.Tecplot;

public class Interpolation {
	public static double[] newtonInterpolationCoef(double[] x, double[] y) {
		int n = x.length;
		double[][] v = new double[n][];
		for(int i = 0; i < n; i++) {
			v[i] = new double[n - i];
		}
		
		for(int i = 0; i < n; i++) {
			v[i][0] = y[i];
		}
		for(int i = 1; i < n; i++) {
			for(int j = 0; j < n - i; j++) {
				v[j][i] = (v[j+1][i-1] - v[j][i-1]) / (x[j + i] - x[j]);
			}
		}
		System.out.println(MVO.toString(v));
		return v[0];
	}
	
	public static double newtonInterpolationValue(double[] x, double[] y, double x0) {
		return MVO.polyval(newtonInterpolationCoef(x, y), x, x0);
	}
	
	
	public static double[][] spline(double[] x, double[] y, int type, double[] param) {
		int n = x.length;
		double[][] coef = new double[n - 1][3];
		TridiagMatrix triMat = new TridiagMatrix(n);
		double[] c = new double[n];
		double[] b = new double[n];
		
		double deltaX0, deltaX1, deltaY0, deltaY1;
		for(int i = 1; i < n - 1; i++) {
			deltaX0 = x[i] - x[i - 1];
			deltaX1 = x[i + 1] - x[i];
			triMat.lower[i - 1] = deltaX0;
			triMat.diag[i] = 2 * (deltaX0 + deltaX1);
			triMat.upper[i] = deltaX1;
			deltaY0 = y[i] - y[i - 1];
			deltaY1 = y[i + 1] - y[i];
			b[i] = 3 * (deltaY1 / deltaX1 - deltaY0 / deltaX0);
		}
		
		if (type == 0) {
			triMat.diag[0] = 1;
			triMat.diag[n - 1] = 1;
			b[0] = 0;
			b[n - 1] = 0;
		} else if (type == 1) {
			triMat.diag[0] = 2;
			triMat.diag[n - 1] = 2;
			b[0] = param[0];
			b[n - 1] = param[1];
		} else if (type == 2) {
			deltaX0 = x[1] - x[0];
			deltaY0 = y[1] - y[0];
			triMat.diag[0] = 2 * deltaX0;
			triMat.upper[0] = deltaX0;
			deltaX1 = x[n - 1] - x[n - 2];
			deltaY1 = y[n - 1] - y[n - 2];
			triMat.lower[n - 2] = deltaX1;
			triMat.diag[n - 1] = 2 * deltaX1;
			b[0] = 3 * (deltaY0 / deltaX0 - param[0]);
			b[n - 1] = 3 * (param[1] - deltaY1 / deltaX1);
		} else if (type == 3) {
			triMat.diag[0] = 1;
			triMat.upper[0] = -1;
			triMat.diag[n - 1] = 1;
			triMat.lower[n - 2] = -1;
		} else {
			throw new IllegalArgumentException("Wrong Spline type!");
		}
		
		
		triMat.LUFactorize();
		triMat.LUSolve(b, c);
		
		for(int i = 0; i < n - 1; i++) {
			deltaX0 = x[i + 1] - x[i];
			deltaY0 = y[i + 1] - y[i];
			coef[i][0] = deltaY0 / deltaX0 - deltaX0 * (2 * c[i] + c[i + 1]) / 3;
			coef[i][1] = c[i];
			coef[i][2] = (c[i + 1] - c[i]) / (3 * deltaX0);
		}
		return coef;
	}
	
	public static double splineValue(double[] x, double[] y, double[][] coef, double x0) {
		double y0;
		int idx = 0;
		for(int i = 1, n = x.length; i < n; i++) {
			if (x0 <= x[i]) {
				break;
			}
			idx++;
		}
		
		double dx = x0 - x[idx];
		y0 = coef[idx][2] * dx;
		y0 = (y0 + coef[idx][1]) * dx;
		y0 = (y0 + coef[idx][0]) * dx + y[idx];
		return y0;
	}
	
	public static void main(String[] args) {
//		double[] x = {0, 2, 3};
//		double[] y = {1, 2, 4};
//		double[] c = newtonInterpolationCoef(x, y);
//		System.out.println(Arrays.toString(c));
//		
//		x = new double[4];
//		y = x.clone();
//		for(int i = 0; i < 4; i++) {
//			x[i] = Math.PI * i / 6.0;
//			y[i] = Math.sin(x[i]);
//		}
//		c = newtonInterpolationCoef(x, y);
//		System.out.println(Arrays.toString(c));
		
//		double[] x = {0, 1, 2};
//		double[] y = {3, -2, 1};
//		double[][] coef = spline(x, y, 0, null);
//		System.out.println(MVO.toString(coef));
		
		double[] x = {0, 1, 2, 3, 4, 5};
		double[] y = {3, 1, 4, 1, 2, 0};
		double[][] coef = spline(x, y, 3, null);
		System.out.println(MVO.toString(coef));
		double[] xx = MVO.linspace(0, 5.0, 100);
		double[] yy = new double[xx.length];
		for(int i = 0, n = xx.length; i < n; i++) {
			yy[i] = splineValue(x, y, coef, xx[i]);
		}
		double[][] yy0 = new double[1][];
		yy0[0] = yy;
		Tecplot.LineXY("mySpine.dat", xx, yy0);
	}
}
